//SUMMARY:  This do file replicates all presented analyses in Ballard-Rosa, Rogowski, Scheve & Thor "Inequality, Information, and Income Tax Policy Preferences in Austria and Germany," PSRM
* PLEASE NOTE:  This file generates all output drawn from individual-level analyses on German data.  This file should be supplemented with the "PSRM replication - Austrian results.do" file to generate full output.

clear all
set more off

//Set classpath:
cd "XXX"

use "Germany (indiv) PSRM replication.dta", clear

svyset caseid [pweight = weight]
set scheme s1mono

** Figure 1 (German half):  Menu estimates of preferred tax rates
keep if treatment == 0 //Only want this figure for respondents in control condition

rename BracketVAT Bracket6
reshape long Bracket, i(caseid) j(index)
rename Bracket rate_choice
rename index bracket
isid caseid bracket

epctile rate_choice, p(50) svy over(bracket)
matrix p50 = r(table)[1,1...]'
epctile rate_choice, p(75) svy over(bracket)
matrix p75 = r(table)[1,1...]'
preserve
	drop if bracket == 1
	epctile rate_choice, p(25) svy over(bracket)
	matrix p25 = (0 \ r(table)[1,1...]')

	clear
	svmat p25, names(loq)
	svmat p50, names(med)
	svmat p75, names(upq)
	gen iqr = upq - loq
	gen bracket = _n
	
	gen rate_curr = 99
	replace rate_curr = 0 if bracket == 1
    replace rate_curr = 18.76 if bracket == 2
    replace rate_curr = 29.95 if bracket == 3
    replace rate_curr = 42 if bracket == 4
    replace rate_curr = 45 if bracket == 5
    replace rate_curr = 19 if bracket == 6
    assert rate_curr != 99
	
	rename (loq1 med1 upq1) (loq med upq)
	tempfile percentiles
	save `percentiles'
restore
merge m:1 bracket using `percentiles', nogen norep assert(3)

egen upper = max(min(rate_choice, upq + 1.5 * iqr)), by(bracket)
egen lower = min(max(rate_choice, loq - 1.5 * iqr)), by(bracket)

collapse (first) upper lower med loq upq rate_curr, by(bracket) fast

 twoway (rbar med upq bracket, horizontal blc(gs10) bfc(ltblue) barw(0.35)) ///
    (rbar med loq bracket, horizontal blc(gs10) bfc(ltblue) barw(0.35)) ///
    (rspike upq upper bracket, horizontal lcolor(gs10)) ///
    (rspike loq lower bracket, horizontal lcolor(gs10)) ///
    (rcap upper upper bracket, horizontal lcolor(gs10)) ///
    (rcap lower lower bracket, horizontal lcolor(gs10)) ///
    (scatter bracket rate_curr, ms(d) mc(maroon*0.7) msize(*1.5) legend(order(7 "Current Rate") bmargin(55 0 15 0) region(lstyle(none) color(none)) ring(0) pos(5)) ///
        yla(1 `" "below" "8,131" "'  2 `" "8,131 -" "13,470 " "' 3 `" "13,471 -" "52,882 " "' ///
        4 `" "52,883 -" "250,731 " "' 5 `" "above" "250,731" "' 6 "VAT", ang(h) notick nogrid) ///
        graphregion(fcolor(white) lcolor(white)) xtitle("") ///
        ytitle("Tax Bracket (Euros)") xsc(r(0 100)) xla(0 "0%" 20 "20%" 40 "40%" 60 "60%" 80 "80%" 100 "100%") ///
        xtitle("") b1title(" ", size(* .025)) b2title("Preferred Tax Rate", color(gs10) size(medium)))

graph export "Output/slider_rates_per_bracket_germany_control.eps", replace

//Figure 2: Illustration of Kakwani Index
*Created by hand, not generated from data.

//Reload full dataset
use "Germany (indiv) PSRM replication.dta", clear

//Define covariate lists
global demog "female age highEduc isMarried hasChildren notEmployed monthlyInc"
global otherCov "hardWork leftVoter"

//TABLE 1 (German half):  Kakwani experimental results
reg kakwani_index i.treatment [pweight=weight]
outreg2 using  "Output/kakwani_combined", replace ctitle(Kakwani (Germany)) tex(frag) nocons bdec(3) label

//FIGURE 3:  Earn more/less treatments by income terciles:
*Income terciles
gen incomeTercile = 1
replace incomeTercile = 2 if monthlyIncome > 3 //Corresponds to monthly income > 1133
replace incomeTercile = 3 if monthlyIncome > 6 //Corresponds to monthly income > 2567
reg kakwani_index i.treatment##incomeTercile [pweight=weight] if treatment != 3 & treatment != 4
margins, dydx(treatment) at(incomeTercile = (1 2 3))
marginsplot, xdim(incomeTercile) yline(0) xtitle(Income tercile) ytitle(Effect on Kakwani index) title("Positional treatment effects, by income tercile") legend(order(4 "Earn less" 3 "Earn more"))
graph export "Output/kakwani(earnXincome).png", as(png) replace


//APPENDIX RESULTS:

//TABLE A.2:  Demographic balance table for German sample:
* Note that population level data were entered by hand to the table generated by the code below to create final table

preserve
    tab gender, gen(g)
    gen a18_34 = inrange(age, 18, 34)
    gen a35_54 = inrange(age, 35, 54)
    gen a55p = age>=55
    tab educ_de, gen(e)

    label var g1 "Gender: Male"
    label var g2 "Gender: Female"
    label var a18_34 "Age: 18–34"
    label var a35_54 "Age: 35–54"
    label var a55p "Age: 55+"
    label var e1 "No secondary general school-leaving certificate"
    label var e2 "Presently attending schools"
    label var e3 "Secondary general school-leaving certificate"
    label var e4 "Certificate of ten-grade school of general education"
    label var e5 "Intermediate school-leaving certificate"
    label var e6 "Fachhochschule/University entrance qualification"
    label var e7 "Bachelor"
    label var e8 "Master"
	label var e9 "Diploma"
	label var e10 "Doctor's degree"

    local rows g1 g2 a18_34 a35_54 a55p e1-e10
	
	foreach v of varlist `rows' {
        replace `v' = 100 * `v'
    }
	
    eststo clear
    estpost summarize `rows' [aw=weight]
    eststo Weighted
    estpost summarize `rows'
    eststo Raw

	esttab Weighted Raw using "Output/sample_statistics_de.tex", ///
       cells("mean(fmt(%3.2f))") mtitle("Weighted Sample" "Raw Sample") label ///
	   replace noobs
restore

//Table A.3:  Manipulation check for positional inequality treatments
* Info on income distribution only provided for Ineq treatment, so we combine earnings groups with control:
gen ineqInfo = 0
replace ineqInfo = 1 if treatment == 3
replace ineqInfo = 2 if treatment == 4
lab def ineqInfoLab 0 "Control/Earnings" 1 "Ineq. 99/50" 2 "Ineq. 99/10"
lab val ineqInfo ineqInfoLab
//Generate "close" estimate of correct answer if respondent was within one category
gen close1pc = (onePCEstimate == 5 | onePCEstimate == 6)
gen closeMedian = (medianEstimate == 3 | medianEstimate == 4 | medianEstimate == 5)
gen closeBottom = (bottom10PCEstimate == 3 | bottom10PCEstimate == 4 | bottom10PCEstimate == 5)

reg close1pc i.ineqInfo [pweight=weight]
outreg2 using  "Output/manip_Germany", replace ctitle(99th pct.) tex(frag) nocons bdec(3) label
reg close1pc i.ineqInfo $demog i.ptTime5 [pweight=weight]
outreg2 using  "Output/manip_Germany", append ctitle(99th pct.) tex(frag) nocons bdec(3) label
reg closeMedian i.ineqInfo [pweight=weight]
outreg2 using  "Output/manip_Germany", append ctitle(Median) tex(frag) nocons bdec(3) label
reg closeMedian i.ineqInfo $demog i.ptTime5 [pweight=weight]
outreg2 using  "Output/manip_Germany", append ctitle(Median) tex(frag) nocons bdec(3) label
reg closeBottom i.ineqInfo [pweight=weight]
outreg2 using  "Output/manip_Germany", append ctitle(10th pct.) tex(frag) nocons bdec(3) label
reg closeBottom i.ineqInfo $demog i.ptTime5 [pweight=weight]
outreg2 using  "Output/manip_Germany", append ctitle(10th pct.) tex(frag) nocons bdec(3) label

//Table A.5:  Manipulation check for earnings treatments
* Info on income percentiles only provided for Earn treatment, so we combine earnings groups with control:
gen earnInfo = 0
replace earnInfo = 1 if treatment == 1
replace earnInfo = 2 if treatment == 2
lab def earnInfoLab 0 "Control/Ineq." 1 "Earn more" 2 "Earn less"
lab val earnInfo earnInfoLab

gen ownIncMore = 100 - ownIncPct
gen incPctDiff = ownIncMore - earnMoreEstimate
gen absPctDiff = abs(incPctDiff)
reg absPctDiff i.earnInfo [pweight=weight]
outreg2 using  "Output/earnManip_Germany", replace ctitle(Own income pct.) tex(frag) nocons bdec(3) label
reg absPctDiff i.earnInfo $demog i.ptTime5 [pweight=weight]
outreg2 using  "Output/earnManip_Germany", append ctitle(Own income pct.) tex(frag) nocons bdec(3) label

//Table A.7:  Kakwani index experimental results, with controls (German half)
reg kakwani_index i.treatment $demog i.time5 [pweight=weight]
outreg2 using  "Output/kakwani_combined_controls", replace ctitle(Kakwani (Germany)) tex(frag) nocons bdec(3) label

//Figure A.1:  Balance table
* In order to make continuous variables comparable, we standardize them:
egen stdAge = std(age)
egen stdInc = std(monthlyInc)

global demogStd "female stdAge highEduc isMarried hasChildren notEmployed stdInc"

foreach var in $demogStd {
	reg `var' i.treatment
	est sto balGer`var'
}
coefplot (balGerfemale, label(Female)) (balGerstdAge, label(Age)) (balGerhighEduc, label(Higher educ.)) (balGerisMarried, label(Married)) (balGerhasChildren, label(Children)) (balGernotEmployed, label(Unemployed)) (balGerstdInc, label(Income)), drop(_cons) xline(0)
graph export "Output/balanceGermany.png", as(png) name("Graph") replace

//Table A.8:  Addressing non-compliance through 2SLS (German half)
gen earnMore = (treatment == 1)
gen earnLess = (treatment == 2)
gen inc9950 = (treatment == 3)
gen inc9910 = (treatment == 4)
ivreg2 kakwani (complier* = earnMore earnLess inc9950 inc9910) [pweight=weight]
outreg2 using  "Output/kakwaniIV_combined", replace ctitle(Kakwani IV (Germany)) tex(frag) nocons bdec(3) label adds(Wald F-stat, e(widstat))

//Table A.9:  Government trust results (German half)
reg trustGovt i.treatment [pweight=weight]
outreg2 using  "Output/govtTrust", replace ctitle(Trust govt. (Germany)) tex(frag) nocons bdec(3) label
gen doTrust = (trustGovt < 3)
reg kakwani_index i.treatment##i.doTrust [pweight=weight]
outreg2 using  "Output/govtTrust", append ctitle(Kakwani (Germany)) tex(frag) nocons bdec(3) label

//Table A.10:  Revenue raised results (German half)
reg revenue_perc_change_sliders i.treatment  [pweight=weight]
outreg2 using  "Output/revenue", replace ctitle(Revenue raised (Germany)) tex(frag) nocons bdec(3) label
reg revenue_perc_change_sliders i.treatment $demog i.time5 [pweight=weight]
outreg2 using  "Output/revenue", append ctitle(Revenue raised (Germany)) tex(frag) nocons bdec(3) label

//CONJOINT RESULTS
use "Germany CJ.dta", clear
svyset caseid [pweight=weight]

//Figure A.3:  Pooled conjoint results (German half)
reg chosePlan i.rate1 i.rate2 i.rate3 i.rate4 i.rate5 i.rate6 i.revenue [pweight=weight], cluster(caseid)
estimate store baseline

coefplot baseline, mcolor(black) ciopts(lcolor(black) lwidth(thin)) drop(_cons) omitted base xline(0) headings(1.rate1 = "{bf:< 8,131}" 1.rate2 = "{bf:8,131-13,470}" 1.rate3 = "{bf:13,471-52,882}" 1.rate4 = "{bf:52,883-250,731}" 1.rate5 = "{bf:> 250,731}" 1.rate6 = "{bf:VAT}" 1.revenue = "{bf:Revenue}")  ///
	ylabel(, labsize(medlarge)) xtitle("Change in Pr(Tax Plan Selected)") title("Pooled sample") xsize(5) ysize(7) scale(.6)
graph export "Output/CJ pooled (Germany).png", as(png) name("Graph") replace

//Figure A.4:  Conjoint results, control group only from framing experiment (German half)
reg chosePlan i.rate1 i.rate2 i.rate3 i.rate4 i.rate5 i.rate6 i.revenue if treatment == 5 [pweight=weight], cluster(caseid)
estimate store baselineControl

coefplot baselineControl, mcolor(black) ciopts(lcolor(black) lwidth(thin)) drop(_cons) omitted base xline(0) headings(1.rate1 = "{bf:< 8,131}" 1.rate2 = "{bf:8,131-13,470}" 1.rate3 = "{bf:13,471-52,882}" 1.rate4 = "{bf:52,883-250,731}" 1.rate5 = "{bf:> 250,731}" 1.rate6 = "{bf:VAT}" 1.revenue = "{bf:Revenue}")  ///
	ylabel(, labsize(medlarge)) xtitle("Change in Pr(Tax Plan Selected)") title("Control") xsize(5) ysize(7) scale(.6)
graph export "Output/CJ pooled, control only (Germany).png", as(png) name("Graph") replace

//Figure A.5:  Variation in tax plan support by brackets (German half)
local reg_command = "reg chosePlan i.rate1 i.rate2 i.rate3 i.rate4 i.rate5 i.rate6 i.revenue [pweight=weight], cluster(caseid)"
`reg_command'
local full_sse = e(rss)

foreach var in rate1 rate2 rate3 rate4 rate5 rate6 revenue {
    local reg_command_temp = subinstr("`reg_command'", " i.`var'", "", .)
    `reg_command_temp'
    local `var'_sse = e(rss)
    local `var'_partial_r_sq = (``var'_sse' - `full_sse') / ``var'_sse'
}

clear
set obs 7
gen bracket = ""
gen partial_r_sq = .

local i = 1
foreach var in rate1 rate2 rate3 rate4 rate5 rate6 revenue {
    replace bracket = "`var'" in `i'
    replace partial_r_sq = ``var'_partial_r_sq' in `i'
    local i = `i' + 1
}

assert bracket != ""
assert partial_r_sq != .

graph bar partial_r_sq, over(bracket, relabel(1 `" "Bottom" "rate" "' 2 "Rate 2" ///
    3 "Rate 3" 4 "Rate 4" 5 `" "Top" "rate" "' 6 "VAT" 7 "Revenue")) ///
    ytitle("Partial R-Squared") graphregion(color(white)) yla(, nogrid)
graph export "Output/partial_r_squared_per_bracket_germany.eps", replace

//Figure A.6:  Aggregate inequality conjoint results (German half)
use "Germany CJ.dta", clear
svyset caseid [pweight=weight]

gen aggIneq = (treatment == 3 | treatment == 4)
reg chosePlan i.rate1 i.rate2 i.rate3 i.rate4 i.rate5 i.rate6 i.revenue [pweight=weight] if aggIneq == 1, cluster(caseid)
estimate store aggIneqResult
reg chosePlan i.rate1 i.rate2 i.rate3 i.rate4 i.rate5 i.rate6 i.revenue [pweight=weight] if treatment == 5, cluster(caseid)
estimate store controlResult

coefplot (controlResult, label(Control) msymbol(circle) mcolor(black) msize(large) ciopts(lcolor(black)) lwidth(thin)) (aggIneqResult, label(Aggregate Inequality) msymbol(circle_hollow) mcolor(black) msize(large) ciopts(lcolor(black)) lwidth(thin)) , ///
	mcolor(black) ciopts(lcolor(black) lwidth(thin)) drop(_cons) omitted base xline(0) headings(1.rate1 = "{bf:< 8,131}" 1.rate2 = "{bf:8,131-13,470}" 1.rate3 = "{bf:13,471-52,882}" 1.rate4 = "{bf:52,883-250,731}" 1.rate5 = "{bf:> 250,731}" 1.rate6 = "{bf:VAT}" 1.revenue = "{bf:Revenue}")  ///
	ylabel(, labsize(medlarge)) xtitle("Change in Pr(Tax Plan Selected)") xsize(5) ysize(7) scale(.6)
graph export "Output/Agg Ineq CJ results (Germany).png", as(png) replace

//Figure A.7:  Self-centered inequality conjoint results (German half)
gen selfIneq = (treatment == 1 | treatment == 2)
reg chosePlan i.rate1 i.rate2 i.rate3 i.rate4 i.rate5 i.rate6 i.revenue [pweight=weight] if selfIneq == 1, cluster(caseid)
estimate store selfIneqResult

coefplot (controlResult, label(Control) msymbol(circle) mcolor(black) msize(large) ciopts(lcolor(black)) lwidth(thin)) (selfIneqResult, label(Self-Centered Inequality) msymbol(circle_hollow) mcolor(black) msize(large) ciopts(lcolor(black)) lwidth(thin)) , ///
	mcolor(black) ciopts(lcolor(black) lwidth(thin)) drop(_cons) omitted base xline(0) headings(1.rate1 = "{bf:< 8,131}" 1.rate2 = "{bf:8,131-13,470}" 1.rate3 = "{bf:13,471-52,882}" 1.rate4 = "{bf:52,883-250,731}" 1.rate5 = "{bf:> 250,731}" 1.rate6 = "{bf:VAT}" 1.revenue = "{bf:Revenue}")  ///
	ylabel(, labsize(medlarge)) xtitle("Change in Pr(Tax Plan Selected)") xsize(5) ysize(7) scale(.6)
graph export "Output/Self Ineq CJ results (Germany).png", as(png) replace

//Figure A.8:  All treatment conditions conjoint results (German half)
reg chosePlan i.rate1 i.rate2 i.rate3 i.rate4 i.rate5 i.rate6 i.revenue [pweight=weight] if treatment == 1, cluster(caseid)
estimate store a1Result
reg chosePlan i.rate1 i.rate2 i.rate3 i.rate4 i.rate5 i.rate6 i.revenue [pweight=weight] if treatment == 2, cluster(caseid)
estimate store a2Result
reg chosePlan i.rate1 i.rate2 i.rate3 i.rate4 i.rate5 i.rate6 i.revenue [pweight=weight] if treatment == 3, cluster(caseid)
estimate store b1Result
reg chosePlan i.rate1 i.rate2 i.rate3 i.rate4 i.rate5 i.rate6 i.revenue [pweight=weight] if treatment == 4, cluster(caseid)
estimate store b2Result

coefplot (a1Result, label(Earn less) msymbol(circle_hollow) mcolor(black) msize(large) ciopts(lcolor(black)) lwidth(thin)) (a2Result, label(Earn more) msymbol(triangle_hollow) mcolor(black) msize(large) ciopts(lcolor(black)) lwidth(thin)) (b1Result, label(99/50) msymbol(diamond_hollow) mcolor(black) msize(large) ciopts(lcolor(black)) lwidth(thin)) (b2Result, label(99/10) msymbol(square_hollow) mcolor(black) msize(large) ciopts(lcolor(black)) lwidth(thin)) (controlResult, label(Control) msymbol(circle) mcolor(black) msize(large) ciopts(lcolor(black)) lwidth(thin)), ///
	mcolor(black) ciopts(lcolor(black) lwidth(thin)) drop(_cons) omitted base xline(0) headings(1.rate1 = "{bf:< 8,131}" 1.rate2 = "{bf:8,131-13,470}" 1.rate3 = "{bf:13,471-52,882}" 1.rate4 = "{bf:52,883-250,731}" 1.rate5 = "{bf:> 250,731}" 1.rate6 = "{bf:VAT}" 1.revenue = "{bf:Revenue}")  ///
	ylabel(, labsize(medlarge)) xtitle("Change in Pr(Tax Plan Selected)") xsize(5) ysize(7) scale(.6)
graph export "Output/CJ results, all treatments (Germany).png", as(png) replace
